# Simple facade for NetCDF I/O (read_ncdf.f90, read_bg.f90, read_lsm.f90, read_orog.f90)

module NetCDFIO

using NCDatasets

"""
    read_ncdf(filename, varname, lonname, latname, timename,
              lono, lato, nx, ny, jd, nread)

Port of read_ncdf.f90. Returns (cnt, flx::Array{Float64,3}) with shape (nx, ny, cnt).
`jd` is a julian-day-like value; we assume the time variable is days since 1900-01-01,
and compute `jdi = jd - juldate(19000101,0)` upstream before calling if needed.
"""
function read_ncdf(filename::AbstractString, varname::AbstractString,
                   lonname::AbstractString, latname::AbstractString, timename::AbstractString,
                   lono::Real, lato::Real, nx::Integer, ny::Integer, jdi::Real, nread::Integer)
    isfile(filename) || error("read_ncdf: cannot find $(filename)")
    ds = NCDataset(filename)
    try
        lon = vec(ds[lonname][:])
        lat = vec(ds[latname][:])
        time = vec(ds[timename][:])
        xres = lon[2] - lon[1]
        yres = lat[2] - lat[1]

        # Checks analogous to Fortran
        (lon[1] - xres/2 <= lono) || error("incompatible longitude origin")
        (lat[1] - yres/2 <= lato) || error("incompatible latitude origin")
        (lono >= minimum(lon .- xres/2)) || error("longitude out of range")
        (lono <= maximum(lon .- xres/2)) || error("longitude out of range")
        (lato >= minimum(lat .- yres/2)) || error("latitude out of range")
        (lato <= maximum(lat .- yres/2)) || error("latitude out of range")

        ix = argmin(abs.(lon .- (lono + xres/2)))
        jy = argmin(abs.(lat .- (lato + yres/2)))
        kt = argmin(abs.(time .- jdi))

        nt = length(time)
        cnt = (kt + nread - 1) > nt ? (nt - kt + 1) : Int(nread)

        var = ds[varname]
        # Use start,count slicing like Fortran: (lon,lat,time)
        flx = Array(var[ix:ix+nx-1, jy:jy+ny-1, kt:kt+cnt-1])
        return cnt, flx
    finally
        close(ds)
    end
end

"""
    read_lsm!(dom, filename, varname, lonname, latname)

Reads a land-sea mask on the inversion grid and stores into dom.lsm.
"""
function read_lsm!(dom, filename::AbstractString, varname::AbstractString,
                   lonname::AbstractString, latname::AbstractString)
    isfile(filename) || error("read_lsm!: cannot find $(filename)")
    ds = NCDataset(filename)
    try
        lon = vec(ds[lonname][:])
        lat = vec(ds[latname][:])
        xres = lon[2] - lon[1]
        yres = lat[2] - lat[1]
        ix = argmin(abs.(lon .- (dom.rllx + dom.rdx/2)))
        jy = argmin(abs.(lat .- (dom.rlly + dom.rdy/2)))
        size_ok = (ix > 0) && (jy > 0) && (length(lon) >= ix + dom.nxregrid - 1) && (length(lat) >= jy + dom.nyregrid - 1)
        size_ok || error("read_lsm!: incompatible domain or dimensions")
        lsmv = ds[varname]
        window = lsmv[ix:ix+dom.nxregrid-1, jy:jy+dom.nyregrid-1]
        dom.lsm = Float64.(coalesce.(window, 0.0))
        return dom
    finally
        close(ds)
    end
end

end # module
